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ABSTRACT 

The  computer  has  given  rise  to  a  new  mode  of  scientific  practice, 
and  today  computational  science  stands  beside  theory  and  experiment 
as  a  fundamental  methodology.  The  impact  of  the  computer  revolu- 
tion on  science  can  be  projected  from  current  trends.  The  demands  to 
be  made  on  computing  methodologies  will  be  reviewed.  One  of  the 
demands  is  an  ongoing  need  for  excellence  in  computational  metho- 
dologies. 

Generic  difficulties  encountered  in  meeting  these  challenges  will 
be  discussed.  Recent  work  of  the  authors  and  others  will  be  reviewed 
in  this  context. 
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1.  The  Computer  Revolation. 

The  ages  of  history  are  demarcated  as  much  by  the  rise  and  fall  of  ideas  and 
technology  as  by  the  rise  and  fall  of  empires  and  nations.  In  Jiis  sense  the  com- 
puter age  is  a  turning  point  for  human  events  as  was  the  earlier  industrial  revolution 
based  upon  mechanical  energy.  With  a  narrower  focus,  we  consider  the  impact  of 
the  computer  revolution  on  science.  Turning  points  in  science  are  marked  by  the 
introduction  of  new  ideas  and  new  tools,  and  especially  by  developments  which  are 
both  tools  and  ideas.  With  the  invention  of  calculus,  the  electron  microscope  and 
DNA  sequencing,  one  finds  an  ability  to  formulate,  understand  and  solve  a  range  of 
problems  which  had  previously  been  inaccessible.  Judged  from  this  perspective,  the 
computer  revolution  will  have  an  impact  on  science  at  least  as  large  as  did  these 
three  examples. 

Very  large  parts  of  science  are  described  in  terms  of  mathematical  equations. 
In  most  cases  the  solutions  to  these  equations  cannot  be  obtained  without  recourse 
to  numerical  computation.  In  the  thirty  years  since  the  advent  of  main  frame  com- 
puters, an  impressive  range  of  problems  have  been  solved  with  the  aid  of  comput- 
ers. Running  a  modern  wind  tunnel  to  test  airfoil  designs  costs  150  million  dollars 
per  year.  Supercomputers  can  explore  a  much  larger  range  of  ideas  than  can  actu- 
ally be  tested,  and  expensive  test  facilities  can  be  reserved  for  the  most  promising 
designs.  In  fact  a  broad  reinge  of  two  dimensional  fluid  flow  problems  are  under 
reasonable  scientific  control,  due  primarily  to  progress  on  numerical  computation  on 
these  problems.  The  list  of  such  successes  could  be  extended  to  considerable  length. 
Lx)oking  to  the  future,  there  can  be  little  doubt  that  the  mathematization  of  various 
subfields  of  biology  is  impending.  As  indications  we  can  cite  the  recent  computer 
assisted  determination  of  the  three  dimensional  structure  of  the  polio  and  cold 
viruses  [37,  49J,  computer  simulation  of  the  heart  with  natural  and  artificial  valves 
[46]  and  algorithms  for  the  analysis  of  DNA  sequences  [33], 

Computers  can  also  be  used  as  an  experimental  tool  to  explore  the  unknown, 
Feigenbaum  discovered  universal  order  in  the  behavior  of  chaos  [20,  17]  through 
computer  experiments,  a  development  which  has  stimulated  and  rejuvenated  several 
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branchcs  of  mathematics  jmd  physics. 

Computers  can  simulate  experimentally  difficult  or  unattainable  parameter 
ranges,  such  as  the  conditions  in  the  interior  of  the  sun  or  a  few  microseconds  after 
the  big  bang  origin  of  the  universe.  They  can  also  simulate  undesirable  parameter 
ranges  such  as  occur  in  safety  studies  to  avoid  accidents  in  chemical  or  nuclear  reac- 
tors. There  have  been  three  types  of  limits  proposed  for  the  ultimate  scope  of  the 
computational  method.  However,  it  is  not  clear  how  firm  these  limits  will  turn  out 
to  be,  and  in  the  authors'  jvdgment,  it  is  no  more  possible  to  set  ultimate  limits  on 
the  scope  of  the  computational  methodology  in  science  than  it  is  to  set  limits  on  the 
scope  of  future  progress  in  theory  or  experiment.  It  has  been  proposed  that  funda- 
mental physical  limits  such  as  the  speed  of  light  will  limit  the  speed  of  future  com- 
puters. However,  the  development  of  parallel  computers,  which  do  many  opera- 
tions concurrently,  may  permit  an  end  run  around  this  problem.  A  class  of  com- 
binatorial problems,  known  as  NP-complete  problems,  arc  effectively  outside  of 
practical  computation.  The  extent  to  which  the  important  problems  which  fall  into 
this  class  can  be  restricted  to  subproblems  which  do  not  lie  in  this  class  or  can  be 
otherwise  modified  so  as  to  be  effectively  computable  is  not  known.  Finally,  it  has 
been  proposed  that  problems  with  producing  complex  computer  code  in  debugged 
and  reliable  form  may  provide  an  outer  limit  for  the  use  of  computers  to  solve  cer- 
tain types  of  problems.  Again  there  are  several  strategics  which  may  mitigate  this 
problem.  We  mention  in  particular:  portable  operating  systems  and  software  tools, 
high  level  languages  and  standardized  calling  sequences  and  modular  libraries  to 
allow  interchangeability  of  reliable  software  components  which  can  be  used  and  thus 
tested  through  time  in  a  variety  of  applications. 

There  is  no  shortage  of  current  limits.  For  all  of  its  astonishing  successes, 
scientific  computations  provide  much  less  than  science  and  technology  need.  In 
most  important  scientific  applications,  the  solutions  arc  undercomputed  relative  to 
the  scientific  requirements.  The  gap  between  needs  and  performance  is  often  sub- 
stantial (and  includes  almost  all  three  dimensional  problems).  Narrowing  this  gap 
will  depend  on  progress  in  both  hardware  and  computational  methodologies.    To 
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asscss  the  current  achievements  and  shortfalls  of  specific  computational  methods  is 
like  trying  to  decide  whether  a  glass  of  water  is  half  full  or  half  empty.  It  can 
perhaps  wryly  be  observed  that  those  with  responsibility  for  supporting  and  main- 
taining a  code  tend  to  believe  that  it  is  correct  and  adequate,  those  who  use  it  tend 
to  believe  that  the  code  is  rather  imperfect,  but  can  nevertheless  serve  as  a  guide  to 
the  wise  scientist  or  engineer,  and  those  who  develop  new  methods  tend  to  be 
acutely  aware  of  the  flaws  in  existing  codes. 

2.   The  Method  of  Scientific  Computing. 

What  arc  the  methodologies  which  define  scientifi';  computing  as  a  distinct 
approach  to  science,  complementing  the  traditional  approaches  of  theory  and  experi- 
ment? Scientific  computing-  begins  with  mathematical  modeling,  whereby  essential 
features  of  a  scientific  problem  are  expressed  in  terms  of  mathematical  equations. 
Typically,  important  compromises  and  judgments  are  made  at  this  stage.  In  order 
to  carry  out  later  computational  steps,  unimportant  and  sometimes  important  prob- 
lem features  must  be  suppressed.  These  judgments  lead  to  a  sequence  of  models 
which  capture  different  aspects  of  the  same  problem.  The  different  models  may  be 
integrated  hierarchically  in  that  a  fine  scale  model  may  be  used  to  set  the  parame- 
ters of  a  coarser  model,  or  they  may  be  integrated  only  in  the  final  judgments  of  the 
scientist  or  engineer  using  them. 

After  the  mathematical  equations  have  been  formulated,  they  must  be  cast  in  a 
form  suitable  for  numerical  computation.  This  usually  requires  that  they  be  discrc- 
tized  in  one  way  or  another.  The  key  step  is  the  formulation  of  a  solution  algorithm 
to  solve  the  discretized  equations.  This  step  is  so  important  that  it  often  influences 
the  preceding  steps  of  model  formulation  and  discretization. 

The  next  stage  is  testing  and  validation  of  the  solution  algorithm.  Validation 
may  be  accomplished  by  the  following  methods:  comparison  with  known  analytic 
solutions,  comparison  with  previously  validated  computations,  comprrison  with 
laboratory  experiments  and  internal  consistency  checks.  Such  internal  checks 
include  convergence  under  mesh  refinement,  ansilysis  of  solution  errors  and  analysis 
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of  diagnostic  data.  Truncation  errors  and  convergence  rates  can  be  analyzed 
mathematically.  Occasionally  there  may  be  a  mathematical  convergence  theorem 
for  the  solution  algorithm. 

The  final  stage  in  scientific  computing  is  to  use  a  validated  and  debugged  code 
for  scientific  or  engineering  purposes. 

The  preceding  outline  for  scientific  computing  would  strike  most  practitioners 
as  glib  as  it  omits  all  discussion  of  difficulties.  We  focus  on  the  difficulties  inherent 
in  the  solution  algorithm  stage.  Dealing  with  these  difficulties  defines  the  subject  of 
numerical  analysis. 

Typical  methodological  difficulties,  as  they  manifest  themselves  to  a  user,  show 
the  following  symptoms:  Slow  convergence,  nonconvergence,  numerical  instabilities 
and  numerical  simulation  of  spurious  physical  effects.  Nonconvergence  is  normally 
the  sign  of  an  incorrect  algorithm,  but  in  listing  it  here,  we  have  in  mind  more  sub- 
tle possibilities,  such  as  solution  convergence  in  an  L  norm,  while  solution  values  at 
specific  points  (Z,^  convergence)  may  be  invalid.  The  trap  is  easy  to  fall  into,  since 
the  solution  value  at  some  singular  point  (such  as  a  crack  tip)  may  be  both  the  most 
important  and  the  least  convergent  part  of  the  solution.  The  discretization  process 
is  a  modification  of  the  mathemf  ical  equations  and  consequcntiy  of  the  physics  or 
problem  formulation  which  the  equations  represent.  In  some  cases  the  original  phy- 
sical process  is  unstable  or  only  weakly  stable  to  changes  in  the  equations  or  equa- 
tion parameters.  In  such  cases,  there  is  the  danger  that  numerical  discretization 
errors  may  grossly  change  the  nature  of  the  solution.  Even  for  stable  physics,  these 
errors  may  change  the  solution  quantitatively  to  an  unacceptable  degree. 

The  symptoms  which  a  user  may  observe  as  a  sign  of  poor  performance  typi- 
cally have  their  origin  in  difficulties  inherent  in  the  original  mathem-tical  equations 
and  their  solutions.  In  general  finite  difference  and  finite  elcmrnt  methods  give 
satisfactory  performance  when  computing  smooth  solutions  to  regular  problems  in  a 
space  of  not  too  high  a  dimension.  We  want  to  emphasize  two  somewhat  unrelated 
t)'pes  of  problem  difficulties.   The  first  has  to  do  with  rapid  solution  variation,  Isttc 
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space  or  time  derivatives  and  jump  discontinuities  or  other  solution  singularities. 
Such  features  are  ubiquitous.  They  arise  in  and  give  rise  to  boundary  layers. 
Chemical  reactions  are  a  typical  source  of  multiple  time  scales,  also  known  as  stiff 
systems.  Problems  with  multiple  material  interfaces  or  shock  waves  have  discon- 
tinuous solutions.  To  explain  this  point  more  fully  we  note  that  the  shock  width  in 
gas  dynamics  is  of  the  order  of  a  mean  free  path,  which  is  about  10~'  cm  and  far 
smaller  than  the  macroscopic  dimensions  of  typical  flow  fields.  The  pressure  gra- 
dient within  the  shock  might  then  be  of  the  order  of  10*  or  more  using  dimensional 
units  appropriate  to  the  problem.  Likewise  a  flame  front  width  could  be  10"^  cm 
but  such  numbers  are  highly  problem  dependent.  For  example  in  an  oil  reservoir 
fire  flood,  the  flame  width  could  be  20  feet,  and  still  small  compared  to  macroscopic 
dimensions.  Chemical  reaction  rates  within  a  single  problem  can  easily  vary  by  a 
factor  of  10^''  because  of  their  exponential  dependence  on  activation  energies.  The 
importance  of  discontinuities  and  singularities  has  long  been  recognized  by  the  com- 
munity working  on  hyperbolic  equations.  However,  they  are  just  as  important  for 
elliptic  problems  such  as  steady  state  elasticity  or  multi-fluid  incompressible  flow, 
due  to  crack  tips,  corners  and  material  or  phase  boundaries,  a  fact  which  appears  to 
have  been  given  less  emphasis  than  it  deserves. 

The  second  major  class  of  problem  difficulties  we  discuss  has  to  do  with  the 
occurrence  of  a  large  or  infinite  number  of  essential  degrees  of  freedom.  Problems 
having  these  features  arise  in  the  study  of  turbulence,  statistical  mechanics  and 
equations  of  state,  quantum  field  theory  and  stochastic  partial  differential  equations. 
This  class  of  difficulties  is  the  more  intractable  of  the  two,  and  the  examples  of 
scientific  studies  [18,  21,  32,  34,  42,  44]  which  can  be  cited  fall  far  short  of  the  sys- 
tematic needs  for  computational  solutions  in  these  areas  due  to  slow  convergence. 
Looking  ahead  to  our  proposed  methodology  for  the  first  class  of  problems,  we 
observe  that  the  problems  with  an  infinite  number  of  degrees  of  freedom  are  also 
very  difficult  from  a  theoretical  point  of  view,  and  that  theoretical  ideas  have  not 
yet  provided  the  means  to  devise  numerical  methods  with  the  required  enhanced 
capabilities  or  convergence  properties. 
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3.   The  Grand  Unified  Scheme  (GUS). 

As  a  general  principle,  we  propose  the  maximal  use  of  analytic  knowledge  of 
solution  properties.  This  principle  has  been  successfully  used  in  the  czsc  of  smooth 
solutions,  where  the  known  solution  regularity  is  used  to  devise  higher  order 
methods  which  give  accelerated  convergence.  For  the  case  of  solution  singularities, 
which  give  rise  to  the  first  class  of  difficulties  mentioned  in  the  previous  section,  we 
propose  maximal  use  of  analytic  knowledge  of  those  solution  singularities.  Perhaps 
paradoxically,  the  solution  singularities  can  be  a  source  of  problem  simplification,  in 
that  the  idealized  or  asymptotic  behavior  in  the  neighborhood  of  the  singularity  may 
possess  a  simpler  structure  or  higher  symmetry  than  does  the  full  solution.  Thus 
the  singularity  may  be  an  opportunity  for  theoretical  progress;  there  may  also  be  a 
considerable  body  of  knowledge  concerning  the  singularities.  Our  proposal  is  to  use 
this  body  of  knowledge  within  the  numerical  algorithm  itself,  and  where  necessary 
to  create  the  required  knowledge.  An  instance  of  this  circle  of  ideas  is  the  vortex 
method,  which  uses  vortices  to  represent  idealized  fluid  elements  in  (for  example) 
shear  layers  and  turbulent  flow  fields  [14,  15,  32,  34,  50,  61,  62].  Related  ideas  arc 
the  interface  methods  based  on  boundary  integral  methods  [3,  58],  and  conformal 
mapping  [43]. 

The  Grand  Unified  Scheme  is  a  version  of  the  general  principles  discussed 
above,  suitable  for  compressible  reactive  fluid  flow.  These  ideas  are  an  outgrowth 
of  work  carried  out  over  several  years  with  many  colleagues,  and  were  presented  in 
preliminary  form  in  [29].  The  scheme  could  also  be  called  the  Kitchen  Smk  Scheme 
(KSS)  because  the  main  idea  is  to  use  everything.  There  are  four  main  components 
to  this  scheme  and  we  discuss  each  briefly  in  this  section.  The  four  components  arc: 
interior  schemes,  front  tracking,  automatic  mesh  refinement  and  automatic  asymp- 
totic reactive  chemistry  (or  automatic  mode  selection). 

High  quality  interior  schemes  have  been  built  around  two  ideas:  flux  limitcrs  to 
control  overshoots  and  numerical  oscillations  and  approximate  Riemann  solvers  to 
achieve  upwind  differencing  [9,  16,  36,  57]. 


The  main  idea  of  front  tracking  is  to  introduce,  as  an  independent  computa- 
tional degree  of  freedom,  the  surfaces  of  jump  discontinuity  which  may  occur  within 
a  solution.  This  method  will  be  explained  in  more  detail  in  the  following  section. 

Automatic  mesh  refinement  is  the  idea  of  using  fine  grids  in  regions  where  the 
solution  is  singular  or  rapidly  varying,  coupled  with  course  grids  in  regions  where 
the  solution  has  a  regular  behavior.  This  method  has  given  rise  to  a  large  ongoing 
enterprise,  from  which  we  cite  representative  works  [1,  2,  8,  56J. 

Automated  asymptotics  or  automated  mode  selection  for  reactive  chemistry  is 
the  proposal  to  determine  the  ambient  conditions  from  the  neighbors  of  a  given 
mesh  block  and  on  the  basis  of  these  conditions  to  determine  rate  limiting  reactions 
and  simplified  chemistry,  so  that  the  rapid  reactions  are  replaced  by  frozen  or 
quasi-steady  state  conditions.  A  hierarchy  of  time  scales  would  be  identified,  and  in 
this  hierarchy,  the  rapid  times  are  treated  as  quasi-steady  state,  the  intermediate 
times  would  be  treated  implicitly,  while  the  slowest  times  would  be  treated  explicitly 
and  the  transients  for  these  time  scales  would  be  fully  and  accurately  resolved. 

There  have  been  several  attempts  to  date  to  implement  this  unified  proposal. 
In  [12]  the  modem  interior  schemes  were  successfully  combined  with  a  limited  ver- 
sion of  front  tracking.  In  [6],  automatic  mesh  refinement  was  successfully  added  to 
the  above.  However,  efforts  to  combine  mesh  refinement  with  modern  interior 
schemes  in  the  absence  of  tracking  led  to  spurious  waves  generated  by  strong  shocks 
at  the  boundaries  of  the  fine  meshes.  For  a  problem  which  does  not  contain  strong 
shocks,  the  successful  combination  of  mesh  refinement  with  a  high  quality  interior 
scheme  has  been  obtained  [7].  The  automated  analysis  of  critical  modes  and  solu- 
tion asymptotics  as  proposed  above  is  an  ambitious  undertaking  and  has  not  been 
attempted  to  date. 

Hybrid  schemes  have  obvious  software  complexity  problems  which  will  be 
mentioned  below.  More  serious  arc  the  scientific  difficiilties.  Here  we  would  like 
to  stress  the  problems  of  a  fundamental  mathematical  nature.  Riemann  problems 
define  the  scale  invariant  large  time  solution  asymptotics  to  leading  order  and  pose  a 


cumber  of  problems  both  for  their  behavior  in  the  large  and  their  behavior  in  higher 
dimensions  [28,  30].  Such  considerations  have  inspired  a  study  of  Ricmann  prob- 
lems in  a  variety  of  physical  contexts  and  have  led  to  the  discovery  of  new  and  strik- 
ing mathematical  phenomena  [19,  51,  52,  53,  54],  including  shocks  which  violate 
most  proposed  physical  admissibility  conditions.  A  related  problem  is  the  hyper- 
bolic problem  with  an  embedded  elliptic  region,  a  notoriously  difficult  problem 
which  arises  in  a  variety  of  contexts  including  transonic  flow,  elasticity,  oil  reser- 
voirs, multi-phase  flow  and  the  relativistic  hydrodynamics  of  a  quark-gluon  plasma 
[4,  5,  38,  55,  59,  60]. 

4,   Front  Tracking 

Front  tracking  can  be  most  adequately  described  by  a  picture.  In  Figure  4.1  wc 
show  the  grids  used  for  a  front  tracking  solution.  Note  that  there  is  a  regular  two 
dimensional  grid  covering  a  channel,  including  the  wedge,  which  is  an  obstacle  in 
the  channel.  The  front  includes  a  shock  wave  incident  on  the  wedge,  the  channel 
walls  and  the  wedge,  as  well  as  the  inlet  and  outlet  surfaces.  It  is  represented  by  a 
lower  dimensional  (one  dimensional)  grid,  which  moves  in  time,  using  velocities 
derived  from  the  jump  conditions  (the  Paemann  problem)  of  gas  dynamics. 

The  plan  to  use  front  tracking  was  initiated  by  O.  McBryan  and  one  of  the 
authors  seven  yc2u-s  ago;  a  clear  and  early  statement  of  the  method  and  the  scientific 
program  is  contained  in  [22].  This  program  was  developed  and  implemented  in  a 
scries  of  more  than  40  papers  by  the  authors,  Oliver  McBryan  and  other  coworkers, 
see  for  example  [10,  11,  13,  22,  23,  24,  25,  26].  In  these  papers,  proof  of  principle 
and  scientific  validation  has  been  established.  Validation  has  been  accomplished  by 
comparison  to  analytic  solutions,  laboratory  experiments,  elementary  one  dimen- 
sional calculations,  and  previously  validated  computer  codes,  as  well  as  convergence 
under  mesh  refinement.  These  tests  have  been  conducted  in  the  context  of  a  variety 
of  different  applications.  In  Figure  4.2  we  show  convergence  under  mesh  refinement 
and  comparison  of  a  two  dimensional  computation  to  a  simpler  one  dimensional 
(radially  symmetric)  computation  for  a  cylindricaUy  symmetrical  detonation  front,  in 
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work  by  B.  Bukict  [10]. 

At  the  outset  the  concept  of  front  tracking  was  considered  to  be  hopelessly 
wrong,  and  the  plan  to  initiate  this  approach  in  a  serious  way  prompted  considerable 
controversy.  The  concerns  were  of  three  types:  This  method  might  not  work  scien- 
tifically, it  was  too  complicated  to  implement,  and  it  might  never  be  able  to  handle 
problems  of  even  moderate  engineering  complexity.  Actually  there  was  a  record  of 
previous  attempts  with  this  method,  which  provided  encouragement  that  the  method 
was  fundamentally  sound  [40,  45,  48].  The  scientific  validation  was  discussed 
above.  Various  numerical  issues  such  as  the  coupling  scheme  between  the  interior 
and  the  front  have  been  examined  but  deserve  further  study.  Local  behavior  in  the 
neighborhood  of  co-dimension  two  intersection  points  of  tracked  discontinuities  also 
needs  further  study. 

The  software  complexity  issues  were  handled  through  a  strategy  which  is  dif- 
ferent from  that  normally  employed  in  computational  fluid  dynamics,  and  which 
might  be  worthy  of  consideration  by  others.  These  issues  have  been  [27]  and  will  be 
discussed  elsewhere. 

On  the  basis  of  this  success  with  scientific  validation  and  control  of  software 
complexity,  it  can  be  stated  that  a  proof  of  scientific  principle  for  front  tracking  has 
been  achieved.  Next  we  discuss  the  transition  of  the  front  tracking  method  from 
benchmark  validation  studies  to  problems  of  moderate  engineering  complexity.  A 
central  research  issue  involved  in  this  transition  is  the  ability  to  allow  bifurcations  or 
changes  in  front  topology  due  to  the  interaction  or  crossing  of  tracked  waves.  This 
program  is  partially  complete  and  we  report  on  progress  achieved  by  coworkers. 

In  Figure  4.3  we  show  a  computation  of  a  line  drive  well  configuration  taken 
from  an  oil  reservoir  study  in  which  the  jet  pinch  off  instability  leads  to  successive 
front  bifurcation  [26]. 

In  Figure  4.4  we  show  successive  stages  in  the  front  topology  as  a  shock  wave 
hits  a  ramp,  forms  a  regular  reflection  and  finally  a  Mach  stem  [11]. 
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In  Figure  4.5  vvc  show  Kclvin-Helmholtz  roll  up  for  compressible  gas  dynam- 
ics, with  the  tracked  slip  line  exiting  and  reappearing  through  periodic  boundaries 
[11]. 

In  Figure  4.6  we  show  successive  stages  in  front  topology  as  a  shock  wave  hits 
a  contact  in  initiation  of  Meshkov  instability,  taken  from  unpublished  work  of  J. 
Grove  [35].  There  arc  two  remarkable  points  in  this  computation.  One  ij  the 
change  in  front  topology  which  ccurs  within  the  computation,  a  transition  which 
has  always  been  regarded  as  a  significant  obstacle  to  the  use  of  this  method.  The 
other  concerns  contiguous  waves,  as  discussed  in  the  next  paragraph. 

There  have  been  several  definitive  achievements  of  the  front  tracking  effort. 
First  wc  mention  a  calculation  of  Grove  [35]  in  which  a  contact  separates  highly 
compressible  (-y  =  1.1)  SF(,  from  air.  When  hit  from  a  shock  on  the  air  side,  the 
transmitted  shock  in  the  SF(,  and  the  contact  itself  lie  almost  on  top  of  each  other, 
making  this  problem  virtually  impossible  to  compute  correctly  by  standard  methods. 
Grove's  solution,  calculated  easily  on  a  20x20  grid  shows  the  efficacy  of  front  track- 
ing. 

J.  Jones  [39]  determined  the  leading  order  endothcrmic  effect  of  radial  cooling 
on  reactive  chemistry  in  curved  detonation  fronts.  This  piece  of  asymptotic  analysis 
solved  a  major  open  problem  in  reactive  hydrodynamics;  it  was  prompted  by  the 
requirement  of  front  tracking  for  analytic  knowledge  of  solution  behavior. 

King,  Lindquist  and  Reyna  [41]  solved  a  20  year  old  problem  posed  by  Ratch- 
ford  [47],  a  leading  petroleum  engineer.  They  showed  that  finite  difference 
methods  can  duplicate  the  phenomena  of  viscous  fingering  with  immiscible  displace- 
ment and  capillary  diffusion,  as  is  observed  in  laboratory  experiments.  This 
achievement  resulted  from  the  need  for  a  benchmark  calculation  for  comparison  to 
front  tracking  studies  [24]. 
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5.  Conclaslons. 

Science  requires  computational  methods  with  enhanced  capabilities  which  go 
beyond  currently  available  techniques.  A  promising  strategy  for  meeting  these 
needs  is  to  use  methods  which  are  adapted  to  known  solution  behavior.  An  outline 
for  pursuing  this  strategy  in  the  context  of  compressible  reactive  flow  has  been 
presented.  Work  directed  at  implementing  this  strategy  has  been  discussed  with  an 
emphasis  on  the  front  tracking  program. 
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7.   Figure  Captions. 

Figure  4.1.  The  grids  used  by  front  tracking  for  a  shock  on  ramp  problem.  A  regu- 
lar two  dimensional  grid  is  superimposed  on  a  one  dimensional  front  (grid).  The 
latter  consists  of  the  incident  and  reflected  shock,  the  ramp  boundary  and  the  boun- 
dary of  the  computational  rectangle. 

Figure  4,2. a.  Comparison  of  one  and  two  dimensional  computations  for  a  cylindri- 
cally  symmetrical  detonation  front  computation,  A  plot  of  pressure  vs.  radius  is 
shown.  The  solid  curve  shows  the  results  obtained  by  the  one-dimensional  random 
choice  computation.  The  vertical  lines  represent  the  range  of  pressure  values  in  the 
two-dimensional  front  tracking  solution  at  a  fixed  radius  as  the  angle  varies.  Thus 
the  vertical  lines  show  the  range  of  angular  dependence  on  the  solution.  The  grid  is 
40  by  40. 

Figure  4.2.b  Convergence  under  mesh  refinement  in  the  two  dimensional  computa- 
tion. 

Figure  4.3.  Automated  bifurcation  of  tracked  fronts  illustrated  in  a  jet  pinch-off  ins- 
tability in  an  oil  reservoir  study.  Here  water  has  been  injected  in  five  wells  located 
at  the  bottom  of  the  figure,  and  oil  is  produced  at  five  wells  located  at  the  top  of  the 
figure.  The  reservoir  parameters  represent  stable  (■.■^placement,  and  so  the  water-oil 
interface  does  not  finger.    However  for  these  same  parameters,  the  oil  is  more 
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mobile,  and  consequently  forms  a  jet  between  the  injection  wells. 

Figure  4.4.  Bifurcations  of  front  topology  in  a  shock  on  ramp  problem.  An  incident 
shock  hits  a  ramp,  and  undergoes  bifurcation  to  a  regular  reflection.  When  the  regu- 
lar reflection  node  reaches  the  top  of  the  ramp,  a  bifurcation  to  a  Mach  type  node 
occurs.  This  Mach  triple  point  is  degenerate  in  that  the  reflected  wave  has  zero 
strength  at  the  triple  point.  The  grid  here  is  30  by  30. 

Figure  4.5.  Compressible  Kelvin-Helmholtz  roll  up,  including  passage  of  the  front 
through  a  periodic  boundary. 

Figure  4.6.  Initiation  of  Mcshkov  instability,  taken  from  unpublished  work  of  J. 
Grove.  In  these  four  frames,  a  shock  wave  hits  a  sinusoidal  contact  discontinuity 
separating  two  gases  of  different  density  (air  and  SF^).  The  first  frame  shows  the 
wave  fronts  just  before  initial  contact,  the  second  after  the  shock  wave  has  hit  the 
contact,  giving  rise  to  transmitted  and  reflected  shocks  in  the  two  gases.  The  third 
and  fourth  figures  show  only  the  main  waves,  consisting  of  the  contact  and  reflected 
and  transmitted  shocks,  after  these  waves  have  separated  from  one  another.  Contin- 
ued to  later  time,  this  nm  will  show  an  interface  instability  known  as  the  Meshkov 
instability. 

Figure  4.7.  The  shock  contact  interaction,  with  highly  compressible  fluid  below  the 
contact,  taken  from  unpublished  work  of  J.  Grove. 
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Fig  4.2,b.      Converge  under  mesh  refinement   in   the   two  dimensional   computation. 
Convergence  of  the  front  and  interior  schemes.  The  pressure  errors  in  the  interior  and 
at  the  front  are  shown  for  NzN  grids  at  the  time  indicated  by  Fig.^%  2b  The  #  signs 
represent  the  interior  error,  where 


Interior  Error  =  100%  x 


lflf\P2.-^Jdxdy 


'0   •'0 
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The  front  error  (error  bars)  gives  the  range  of  the  errors  at  the  front,  defined  as 


Front  Error  =  100%  x 


[P] 


where  [P]  is  the  pressure  jump  at  the  front  in  the  one-dimensional  computation  at  the 
same  time.   The  asterisks  represent  the  error  of  the  average  pressure  behind  the  front, 
lamely 


P-  —P. 

Front  Error(average  pressure)  =  100%  x        '^'^ ~ 
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Figure  4.3 
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Figure   4.4 
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